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Abstract — This paper discusses Hamel's formalism and its 
applications to structure-preserving integration of mechanical 
systems. It utilizes redundant coordinates in order to elim- 
inate multiple charts on the configuration space as well as 
nonphysical artificial singularities induced by local coordinates, 
while keeping the minimal possible degree of redundancy and 
avoiding integration of differential-algebraic equations. 

I. INTRODUCTION 

This paper introduces a new variational integrator for a 
spherical pendulum. The configuration space for this pendu- 
lum is a two-dimensional sphere. Calculations in spherical 
coordinates are not a good option because of unavoidable 
artificial singularities introduced by these coordinates at the 
poles. In addition, the topology of a sphere makes it impos- 
sible to use global singularity-free intrinsic coordinates. 

In order to avoid the issues mentioned above, an integrator 
that utilizes the interpretation of a sphere as a homogeneous 
space was introduced in [15]. This integrator performs very 
well, but has a somewhat large degree of redundancy. This 
paper targets the development of an integrator whose perfor- 
mance is similar to that of the integrator in [15] and whose 
redundancy is the minimum possible. 

Both the present paper and [15] interpret the pendulum 
as a rotating rigid body. The algorithm introduced in [15] is 
based on the evaluation of the rotation matrix that represents 
the attitude for this body. The key feature of the dynamics 
utilized in the present paper is that, in order to capture the 
orientation of a rigid body, it is sufficient to evaluate just 
one column of that rotation matrix. The resulting equations 
of motion are interpreted as Hamel's equations written in 
redundant coordinates. 

The general exposition of discrete Hamel's formalism will 
be a subject of a future publication. Here we demonstrate 
the usefulness of some of this formalism by constructing 
an integrator for a spherical pendulum that is energy- and 
momentum-preserving. The calculations are carried out in 
the Cartesian coordinates of the three-dimensional Euclidean 
space. This allows one to avoid singularities and/or multiple 
coordinate charts that are inevitable for calculations on a 
sphere. Hamel's approach allows one, among other things, 
to represent the dynamics in such a way that the length 
constraint becomes unnecessary. Thus, one avoids the well- 
known difficulty of numerically solving differential-algebraic 
equations. 
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The paper is organized as follows. Hamel's formalism 
and its discretization are briefly discussed in Sections [TT] 
and III The dynamics of a spherical pendulum is reviewed 



in Section IV The discrete model for the pendulum based on 
Hamel's formalism, its comparison to some other discretiza- 
tion techniques, and simulations are given in Sections [V VI 
and|Vlll 

II. LAGRANGIAN MECHANICS 

A. The Euler-Lagrange Equations 

A Lagrangian mechanical system is specified by a smooth 
manifold Q called the configuration space and a function 
L : TQ — > R called the Lagrangian. In many cases, 
the Lagrangian is the kinetic minus potential energy of the 
system, with the kinetic energy defined by a Riemannian 
metric and the potential energy being a smooth function on 
the configuration manifold Q. If necessary, non-conservative 
forces can be introduced (e.g., gyroscopic forces that are 
represented by terms in L that are linear in the velocity), but 
this is not discussed in detail in this paper. 

In local coordinates q = (q 1 , . . . , q n ) on the configuration 
space Q we write L = L(q,q). The dynamics is given by 
the Euler-Lagrange equations 

d_dL = dL . = i 
dt dq l dq l ' 

These equations were originally derived by Lagrange [14] 
in 1788 by requiring that simple force balance F = ma 
be covariant, i.e. expressible in arbitrary generalized coordi- 
nates. A variational derivation of the Euler-Lagrange equa- 
tions, namely Hamilton's principle (also called the principle 
of critical action), came later in the work of Hamilton [12] 
and [13] in 1834/35. For more details, see [4], [18], and 
Theorem 12.11 below. 



B. The Hamel Equations 

In this paragraph we briefly discuss the Hamel equations. 
The exposition follows paper [5]. 

In many cases the Lagrangian and the equations of motion 
have a simpler structure when written using velocity compo- 
nents measured against a frame that is unrelated to system's 
local configuration coordinates. An example of such a system 
is the rigid body. 

Let q = (q 1 , . . . , q n ) be local coordinates on the config- 
uration space Q and u.i £ TQ, i = 1 , . . . , n, be smooth 
independent local vector fields defined in the same coordi- 
nate neighborhood^ The components of Ui relative to the 

'in certain cases, some or all of can be chosen to be global vector 
fields on Q. 



basis d/dq> will be denoted ijjj; that is, 



Ui(q) = i>i(q) 



d_ 



where i,j = 1, . 
stood. 

Let £ = (£\ 
velocity vector q g TQ relative to the basis U\ 

<i = Cui(q); 

then 

l(q,0 := L(q,C Ui (q)) 



, n and where summation on j is under- 
. , £ n ) £ M. n be the components of the 

, U n , I.e. , 

(2) 



(3) 



is the Lagrangian of the system written in the local coor- 
dinates (g, £) on the tangent bundle TQ. The coordinates 
(q, £) are Lagrangian analogues of non-canonical variables 
in Hamiltonian dynamics. 

Define the quantities c™(q) by the equations 



[ui(q),Uj(q)} 



){q)u m {q), 



(4) 



i,j,m = l,...,n. These quantities vanish if and only if 
the vector fields Ui(q), i = 1, ...,n, commute. Here and 
elsewhere, [ • , • ] is the Jacobi-Lie bracket of vector fields 
on Q. 

Viewing it, as vector fields on TQ whose fiber components 
equal (that is, taking the vertical lift of these vector fields), 
one defines the directional derivatives Ui[l] for a function 
I : TQ -> K by the formula 

The evolution of the variables (q, £) is governed by the 
Hamel equations 

d dl m dl 



-e+ Uj [i]. 



(5) 

equations Q 



coupled with equations If it, = d/dq 
become the Euler-Lagrange equations ([TJ. 

Equations |5]) were introduced in [10] (see also [21] for 
details and some history). 

C. Hamilton 's Principle for Hamel 's Equations 

Let 7 : [a, b] — > Q be a smooth curve in the configuration 
space. A variation of the curve j(t) is a smooth map f3 : 
[a,b] x [—s,e] — > Q that satisfies the condition /3(t,0) — 
7(i). This variation defines the vector field 



hit) 



ds 



along the curve j(t). 

Theorem 2.1: Let L : TQ — > K be a Lagrangian and 
I : TQ — > K be its representation in local coordinates (q, £). 
Then, the following statements are equivalent: 
(i) The curve q{t), where a < t < b, is a critical point of 
the action functional 

rb 

L{q, q) dt (6) 



on t\ie space of curves in Q connecting q a to q^ on the 
interval [a, b], where we choose variations of the curve 
q(t) that satisfy 5q(a) — 5q(b) = 0. 
(ii) The curve q(t) satisfies the Euler-Lagrange equa- 
tions ([TJ. 

(hi) The curve (q(t), is a critical point of the functional 

i-b 

l(q,0dt (7) 



with respect to variations St;, induced by the variations 
5q = rfiii{q), and given by 

5? = fi k + 4(g)CV (8) 

The curve (q(t),£(t)) satisfies the Hamel equations |5]) 
coupled with the equations q — £, l Ui{q). 



(iv) 



For the proof of Theorem 2.1 and the early development 



and history of these equations, as well as other variational 
structures associated with Hamel's equations see [22], [10], 
[5], and [2]. 

III. THE DISCRETE HAMEL EQUATIONS 

A. Discrete Hamilton's Principle 

A discrete analogue of Lagrangian mechanics can be 
obtained by discretizing Hamilton's principle; this approach 
underlies the construction of variational integrators. See 
Marsden and West [19], and references therein, for a more 
detailed discussion of discrete mechanics. 

A key notion is that of the discrete Lagrangian, which is 
a map L d : Q x Q —> M that approximates the action integral 
along an exact solution of the Euler-Lagrange equations 
joining the configurations qk 7 qk+i € Q, 

L d (q k ,q k+1 )^ ext / L(q,q)dt, (9) 
q ec([Q,h],Q) J Q 

where C([0, h], Q) is the space of curves q : [0, h] Q with 
q(0) = qk, q(h) — qu+i, and ext denotes extremum. 

In the discrete setting, the action integral of Lagrangian 
mechanics is replaced by an action sum 



S d (q ,qi, ...,9jv) 



N-l 

E 

fc=0 



L d {qk, <7fc+i) ; 



where q k € Q, k — 0, 1, . . . , N, is a finite sequence 
in the configuration space. The equations are obtained by 
the discrete Hamilton's principle, which extremizes the 
discrete action given fixed endpoints go an d Qn- Taking 
the extremum over qi, . . . , qw-i gives the discrete Euler- 
Lagrange equations 

D 1 L d (q k ,q k+1 ) + D 2 L d {q k ^ 1 ,q k ) = 0, 

for k = 1, . . . , N— 1. This implicitly defines the update map 
$ : Q x Q ^ Q x Q, where $(9k-i,?k) = (q k ,q k+ i) and 
Q x Q replaces the velocity phase space TQ of Lagrangian 
mechanics. 

In the case that Q is a vector space, it may be convenient 

to use (q k+ i/ 2 ,Vk,k+i), where q k+l/2 = \{qk + qk+i) and 
Vk,k+i — jiiqk+i — qk), as a state of a discrete mechanical 



system. In such a representation, the discrete Euler-Lagrange 
equations become 

5 (DiL d (q k _ 1/ 2, Vk-i,k) + D 1 L d (q k+1/2 , v k ,k+l)) 
+ } l (D 2 L d (q k _ 1/2 ,v k ^ 1 , k ) - D 2 L d (q k+1/2 ,v k , k+1 )) = 0. 
These equations are equivalent to the variational principle 

N-l 



SS d = {D l L d {q k+1/2 ,v Kk+l )5q^ 



k+1/2 



k=0 



+ D 2 L d (q k+1/2 ,v ktk+1 )8v k , k+1 ) = 0, (10) 

where the variations Sq k+ i/ 2 and 5vk,k+i are induced by the 
variations 5q k and are given by the formulae 

<%+l/2 = \ (<5<7fc+l + Sq k ) , Sv kyk+1 = i (Sq k+1 - Sq k ) . 
B. Discrete Hamel's Equations 

In order to construct the discrete Hamel equations for a 
given mechanical system, one starts by selecting the vector 
fields ui(q), . . . , u n (q) and computing the Lagrangian l(q,£,) 
given by ([3]). One then discretizes this Lagrangian (we only 
discuss the mid-point rule here) and obtains 

' d (?fc+i/27 6c,fc+i) = hl(q k+1 / 2 , £, kik+ i) (11) 

Note that the discretization in ( fTT| is carried out after writing 
the continuous-time Lagrangian as a function of 

One of the challenges of discretizing the Hamel equations 
has been understanding the discrete analogue of the bracket 
term in (j5J. Until recently, it was only known how to handle 
this for systems on Lie groups (see e.g. [6] and [17]). In 
the discrete model of a spherical pendulum discussed below, 
these terms vanish, and we will not discuss the approach 
to discretize the bracket terms in this paper (details on this 
topic can be found in [3]). 

The analogue of the variational principle ( fTO) is obtained 
by setting 

5 ll+i/2 = Wj{<lk+i/i){rf k+l + rf k ), 

m, k+ i = Uvi+i vi) + si 

where B k is the discrete analogue of the bracket term in ([SJ. 
The discrete Hamel equations read 

\{DJ d (qk-i/2iik-i,k) + D u l d (q k+1 / 2 ,£ k , k+ i)) 
+ B* k + ±{D2l d {q k - 1/2 ,Zk-i,k) 

-D 2 l d (q k+1/2 ,^ k+1 ))=0, (12) 

where D u l d is the directional derivative given by the formula 



D u l (qk+i/2, £fc,fc+i) 
d_ 

ds „_, 



+ su(g fe+ i/ 2 ),6,fe+i) 



and where B^ is the discrete analogue of the bracket term 
of the continuous-time Hamel equations Q. As mentioned 
above, this term vanishes for the spherical pendulum prob- 
lem, and is therefore not explicitly shown here (see [3] for 
details). Equations ( fL?) along with the discrete analogue 
of equations (|2]i define the update map (q k -\/ 2 , £fc-i.fc) H> 

{lk+l/2, 6c,fc+l)- 



IV. THE SPHERICAL PENDULUM 

Consider a spherical pendulum whose length is r and mass 
is m. We view the pendulum as a point mass moving on 
the sphere of radius r centered at the origin of M 3 . The 
development here is based on the representation 



M = M x n 

f = r x o 



(13) 
(14) 



of the dynamics of a spherical pendulum, where the pendu- 
lum is viewed as a rigid body rotating about a fixed point. 
Here ft is the angular velocity of the pendulum, M is its 
angular momentum, T is the unit vertical vector (and thus 
the constraint ||r|| = 1 is imposed), and T is the torque 
produced by a force acting on the pendulum, all written 
relative to the body frame. Throughout the paper, boldface 
characters represent three-dimensional vectors. Note that the 
projection of T on T is zero. We assume here that the force is 
conservative, with potential energy U(T). For the pendulum, 
C/(r) = mg(a, T) = mgrT 3 , where a is the vector from the 
origin to the pendulum bobj^Note that the potential energies 
for forces like gravity are invariant with respect to rotations 
about T. 

There are two independent components in equation ( fT3] >. 
We emphasize that this representation, though redundant, 
eliminates the use of local coordinates on the sphere, such 
as spherical coordinates. More details on this appear below. 
Spherical coordinates, while being a nice theoretical tool, 
introduce artificial singularities at the north and south poles. 
That is, the equations of motion written in spherical coor- 
dinates have denominators vanishing at the poles, but this 
has nothing to do with the physics of the problem and is 
solely caused by the geometry of the spherical coordinates. 
Thus, the use of spherical coordinates in calculations is not 
advisable. 

Another important remark is that the length of the vector T 
is a conservation law of equations ( fL3| and ( [14) , 



Til = const, 



(15) 



and thus adding the constraint ||r|| = 1 does not result in 
a system of differential-algebraic equations. The latter are 
known to be a nontrivial object for numerical integration. 

Equations ( fT3j ) and ( fT4] > may be interpreted in a number 
of ways. For instance, one can view them as the dynamics 
of a degenerate rigid body. For this interpretation, select an 
orthonormal body frame with the third vector aligned along 
the direction of the pendulum. The inertia tensor relative to 
such a frame is I — diag{ mr 2 , mr 2 , 0}, and the Lagrangian 
reads 



i(n,r) = Uin,n) - u(t). 



(16) 



With this frame selection, the third component of the angular 
momentum of the body vanishes, 

dl 



Ah 



0. 



2 A11 frames in this paper are orthonormal, and thus the dual vectors are 
interpreted as regular vectors, if necessary. 



and thus there are only two nontrivial equations in ([13). Thus, 
one needs five equations to capture the pendulum dynamics. 
This reflects the fact that rotations about the direction of the 
pendulum have no influence on the pendulum's motion. The 
dynamics then can be simplified by setting SI 3 = 0. 

Alternatively, ( fT3| > and ( fl4| > may be interpreted as the 
dynamics of the Suslov problem (see [21] and [4]) for a 
rigid body with a rotationally-invariant inertia tensor and 
constraint f2 = 0. 

Using either interpretation, the dynamics is represented by 
the system of five first-order differential equations 



or, in components, 



Mi = mgrT 2 , 

f 1 = ~n 2 r 3 , 



AI2 = —mgrT 1 , 

f 2 = r^r 3 , r 3 = ^r 1 



wr 2 



(17) 
(18) 



The latter equations are in fact Hamel's equations written 
in the redundant coordinates (r 1 ,r 2 ,r 3 ) relative to the 
frame 



u 2 



ar 2 
sr 3 



, d 

dr 3 ' 



dv 1 



Recall that the length of T is the conservation law, so that 
the constraint ||r|| = 1 does not need to be imposed, but 
the appropriate level set of the conservation law needs to be 
selected. 

Our discretization will be based on this point of view, i.e., 
the discrete dynamics will be written in the form of discrete 
Hamel's equations. The discrete dynamics will posses the 
discrete version of the conservation law ( fl5] >, so that the 
algorithm should be capable, in theory, of preserving the 
length of r up to machine precision. 



V. VARIATIONAL DISCRETIZATION FOR THE 
SPHERICAL PENDULUM 

The integrator for a spherical pendulum is constructed by 
discretizing Hamel's equations ( fT7| ) and ( p"8j ). 

Let the positive real constant h be the time step. Applying 
the mid-point rule to ( fTo} , the discrete Lagrangian is com- 
puted to be 



2^ ^fc,fc+l) ^k,k+l) 



hU(T k+1/2 ). (19) 



Here f2 fcife+1 = (Sl k k+1 , ft 2 k+1 , 0) is the discrete analogue 
of the angular velocity ft = (fl-,f2 2 ,0) and I\. +1 / 2 = 
!(r&_|_i + F k ). The discrete dynamics then reads 



\ T- (S^fc.fc+i — ^fc-i.fc) 



-fc+l/2> 



(20) 



H 1 'fc+1/2 



^-1/2) — 5 (r^+i/2 + r fe _ 1 / 2 ) 

x i(n M+ i + n fc _i lfc ), (2i) 









-l,k 




(^fc,fc+l 




-l,k 


1 

h 




1 fc- 


1/2 



fc+1/2 



■ k-l/2)i 



(22) 



1 ( V 2 

h V fc+1/2 



ifr 3 
h I 1 fc+1/2 



rfc-1/2) 



-53K r fc+i/2 + rL 1/2 ),(23) 
— l(^fc,fc+l + ^fc-l,fc) 

X (^fc + 1/2 + Tfc-1/2)' 
\ (^l.k+1 + ^fc-l,fc) 

X (^fc+1/2 + Tfc-1/2)' 
1 (^fc,fc+l + ^fc-l,fc) 
X (^fc+1/2 + Tfc-1/2) 
— l(^fc,fc+l + ^fe-l,fe) 
X ( r fc+l/2 + r fc-l/2)' 

(24) 

The discretization of equation ( fl4| ) is constructed to be 
|| r|| -preserving. Indeed, using the isomorphism f2 i-> il 
between the spaces of three-dimensional vectors Q = 
(fl 1 , ri 2 , O 3 ) and skew-symmetric matrices 








-n 3 


n 2 ' 


n = 


n 3 





-n 1 




-VI 2 


fi 1 






equation ( |2"T] i becomes 

r fc+ i/2 = (J - ^r'U + A fe )r fc _ 1/2 , 

where 

^4fc = \ (^fc,fc+i + ^fc-i,fc)|^] 
It is straightforward to check that the matrix 

(i-A^ii + Ak) 

is orthogonal (it is simply the Cayley transform of A k ), and 

therefore ||r fc+ i|| = ||r fc ||. 

As expected, the discrete dynamics is momentum- 
preserving: 

mr2 ( r fc+l/2^L+l + r fc+l/2^fc,fc+l) = C011St ' 

i.e., the vertical component of spatial momentum is con- 
served. This can be verified either using general symmetry 
arguments, or by a straightforward calculation. 
The discrete dynamics is energy-preserving: 



^mr 2 



(^ fc+1 ) 2 + (nl k+1 y 



mgrT 



fc+1/2 



const. 



This is confirmed by multiplying equations (|22) and ( |23] > by 

( fi fc,fc+i+^fc-i,fc) and ( fi ifc+i+ fc-i,fc)' res P ectivel y- and 

adding the result to equation p4| . 

To recap, the proposed method preserves the symplectic 
structure, the length constraint, and the momentum, so by 
a result due to Ge and Marsden [9], the proposed method 
recovers the exact trajectory, up to a possible time reparam- 
eterization. 

3 The matrix I — is invertible if h is sufficiently small. 



VI. COMPARISON WITH OTHER METHODS 

The proposed method takes advantage of the homogeneous 
space structure of S 2 , which has a transitive Lie group action 
by SO (3). In particular, the vector T E S 2 is updated by the 
left action of a rotation matrix, given by the Cayley trans- 
formation of a skew-symmetric matrix Ak that approximates 
the angular momentum ft integrated over a half-timestep. 
Interestingly, this falls out naturally from discretizing the 
Hamel formulation of the spherical pendulum, and it would 
be interesting to see what general choices of coordinate 
frames in the Hamel formulation lead to similar methods for 
more general flows on homogeneous spaces. We will now 
discuss some alternative methods of simulating the spherical 
pendulum equations. 

Homogeneous Space Variational Integrators 

If one were to instead formulate the spherical pendulum 
problem directly on S 2 , it is possible to lift the variational 
principle on S 2 to SO(3), by relating the curve T(t) € 
S 2 with a curve R(t) € SO(3), by the relation T(t) = 
R(t)r(0), where R(0) = I, except that the resulting varia- 
tional principle on SO (3) does not have a unique extremizer, 
due to the presence of a nontrivial isotropy subgroup. With 
a suitable choice of connection, this ambiguity can be 
eliminated, and the resulting problem (and similar problems 
on homogeneous spaces) can be solved using Lie group 
variational integrator techniques, as described in [15]. 

Nonholonomic Integrators 

As mentioned in Section |IV| the spherical pendulum 
equations can be viewed as a Suslov problem, which is 
an example of a nonholonomic mechanical system with no 
shape space. In principle, one could apply a nonholonomic 
integrator, such as the one described in [8] and [20], but 
replacing the length constraint with an infinitesimal con- 
straint and a discrete nonholonomic constraint may result 
in poor numerical preservation of the constraint properties if 
the discrete nonholonomic constraint is poorly chosen. An 
alternative approach to simulating nonholonomic mechanics 
involves a discretization of the forces of constraint, and a 
careful choice of force discretization has been shown to yield 
promising results, see [16] for details. 

Constrained Symplectic Integrators 

Given the relatively simple nature of the unit length 
constraint, it is quite natural to apply the RATTLE algo- 
rithm [1], which is a generalization of the Stormer-Verlet 
method for constrained Hamiltonian systems that is designed 
to explicitly preserve holonomic constraints. This method 
does require the use of a nonlinear solver on a system of 
nonlinear equations of dimension equal to the number of 
constraints. The cost of the nonlinear solve can increase 
significantly as the number of copies of the sphere in the 
configuration space increases. 



Differential-Algebraic Equation Solvers 

The proposed discrete Hamel integrator can be easily 
scaled to an arbitrary number of copies of the sphere, 
possibly chained together in a n-spherical pendulum. Such 
multi-body systems however pose significant challenges for 
differential-algebraic equation solvers, since they are exam- 
ples of what are referred to as high-index DAEs, for which 
the theory and numerical methods are much less developed. 
It is possible to perform index reduction on the system of 
differential-algebraic equations, but this involves significant 
effort, and the numerical results can be mixed. 

Numerical Comparisons 

Since the method is a second-order accurate symplectic 
method, it is natural to compare it to the Stormer-Verlet 
method, as well as the RATTLE method (which is a gen- 
eralization of Stormer-Verlet for constrained Hamiltonian 
systems). 

For the Stormer-Verlet method, we compute Hamilton's 
equations for the spherical pendulum, which is given by 

* = £P. ( 25 ) 
p = -mge 3 + (mgx ■ e 3 - ^||p|| 2 ) x = f(x,p), (26) 

and we apply the generalization of the Stormer-Verlet 
method for general partitioned problems (see (3.4) in [11]), 

Pn+1/2 = Pn + |/(a5n,Pn+l/2), 
X n +1 — X n + — P n +l/2, 
Pn+1 = Pn+1/2 + |/( a; n+l,Pri+l/2)- 

This system of equations is linearly implicit, since the first 
equation is implicit in p„ + i, but the rest of the equations 
are explicit. 

The RATTLE method (see (1.26) in [11]) can be applied 
to the particle in a uniform gravitational field problem, 

x = —p. 

p = -mge 3 , 

subject to the constraint <j>{x) = \(\\x\\ 2 — 1) = 0. We also 
introduce *(cc) = || = x T . Then, the RATTLE method 
applied to this problem is given by, 

Pn+1/2 = Pn - | (mge 3 + *(a;„) T A„) , 

X n+ i = X n + 

= <p(x n+1 ), 
Pn+i = Pn+i/2 - \ {mge 3 + *(x n+1 ) T /i„), 

= i*K+i)'P„+i. 

VII. SIMULATIONS 

In Figures[T]and|2] we present simulations using our theory 
developed above, which we compare with simulations using 
the generalized Stormer-Verlet method and the RATTLE 
method in Figures [3] and |4] respectively. 

For simulations, we select the parameters of the system 
and the initial conditions to be rn = 1 kg, r = 9.8 m, 
h = .2 s, fll = .6 rad/s, fig = rad/s, Tj = .3 m, T§ = 



.932738 m. The trajectory of the bob of the 



.2 m, : 

pendulum is shown in Figure la As expected, it reveals the 
quasiperiodic nature of pendulum's dynamics. Theoretically, 
if one solves the nonlinear equations exactly, and in the 
absence of numerical roundoff error, the Hamel variational 
integrator should exactly preserve the length constraint, and 
the energy. In practice, Figure lb demonstrates that ||r|| stays 
to within unit length to about 10~ 10 after 10,000 iterations. 



Figure lc demonstrates numerical energy conservation, and 
the energy error is to about 10~ 10 after 10,000 iterations 
as well. Indeed, one notices that the energy error tracks the 
length error of the simulation, which is presumably due to 
the relationship between the length of the pendulum and the 
potential energy of the pendulum. The drift in both appear 
to be due to accumulation of numerical roundoff error, and 
could possibly be reduced through the use of compensated 
summation techniques. 




equator. This simulation demonstrates the global nature of 
the algorithm, and also seems to do a good job of hinting at 
the geometric conservation properties of the method. 




Fig. 2: A trajectory with initial conditions above the equator 
integrated with the Hamel integrator. 



We also simulate the spherical pendulum using the gen- 
eralized Stormer-Verlet method and the RATTLE method 



described in Section VI The generalized Stormer-Verlet 
method exhibits surprisingly good unit length preservation in 
Figure 3b of 10~ n when applied to index-reduced version 
of the equations of motion (|25]l-(|26]i. The energy behavior 




(a) Trajectory of the pendulum on S 
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(c) Conservation of energy 

Fig. 1: Numerical properties of the Hamel integrator. 




(a) Trajectory of the pendulum on S 
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Figure [2] shows pendulum's trajectory that crosses the 



Fig. 3: Numerical properties of the Stormer-Verlet method. 



in Figure 3c is typical of a symplectic integrator, with the 
characteristic bounded energy oscillations. Even though the 
RATTLE method is intended to explicitly enforce the unit 
length constraint, it exhibits a unit length preservation in 
Figure Kb] of 1CT 7 , which is poorer than both the Hamel 
variational integrator and the generalized Stormer-Verlet 
method. 



The energy error for RATTLE in Figure 4c 



comparable to that of the generalized Stormer-Verlet method, 
but both pale in comparison to the energy error for the Hamel 
variational integrator. 
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Fig. 4: Numerical properties of the RATTLE method. 

VIII. CONCLUSIONS AND FUTURE WORK 
In this paper, we constructed a variational integrator for the 
spherical pendulum by discretizing Hamel's equations. We 
showed the integrator preserves key mechanical quantities 
and illustrated the work with simulations, and comparisons 
with the generalized Stormer-Verlet method and the RAT- 
TLE method. 

Future work will include the simulation of linked rigid 
body systems as well as the control of such systems. The 
excellent numerical properties of the proposed Hamel varia- 
tional integrator will serve as an excellent basis for construct- 
ing numerical optimal control algorithms, which are heavily 



dependent on the quality of the numerical discretization of 
the natural dynamics. 
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